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Abstract 

A framework for performing dynamic mesh adaptation with the discon- 
tinuous Galerkin method (DGM) is presented. Adaptations include modifi- 
cations of the local mesh step size (/i-adaptation) and the local degree of the 
approximating polynomials (jo-adaptation) as well as their combination. The 
computation of the approximation within locally adapted elements is based 
on projections between finite element spaces (FES), which are shown to pre- 
serve an upper limit of the electromagnetic energy. The formulation supports 
high level hanging nodes and applies precomputation of surface integrals for 
increasing computational efficiency Error and smoothness estimates based 
on interface jumps are presented and applied to the fully /zp-adaptive simu- 
lation of two examples in one-dimensional space. A full wave simulation of 
electromagnetic scattering form a radar reflector demonstrates the applica- 
bility to large scale problems in three-dimensional space. 
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1. Introduction 

The discontinuous Galerkin method [H, 0] nowadays is a well-established 
method for solving partial differential equations, especially for time-dependent 
problems. It has been thoroughly investigated by Cockburn and Shu as well 
as Hesthaven and Warburton, who summarized many of their findings in [i[ 
and respectively. Concerning Maxwell's equations in time-domain, the 
DGM has been studied in particular in The former two apply tetra- 

hedral meshes, which provide flexibility for the generation of meshes also 
for complicated structures. The latter two make use of hexahedral meshes, 
which allow for a computationally more efficient implementation j^. 

In ji[ the authors state that the method can easily deal with meshes with 
hanging nodes since no inter-element continuity is required, which renders it 
particularly well suited for /ip-adaptivity. Indeed, many works are concerned 
with h-, p- or /ip-adaptivity within the DG framework. The first published 



work of this kind is presumably [10| , where the authors consider linear scalar 



hyperbolic conservation laws in two space dimensions. For a selection of 



other publications see |11H15| and references therein. The latter three are 
concerned with the adaptive solution of Maxwell's equations in the time- 
harmonic case. 

In this article, we are concerned with solving the Maxwell equations for 
electromagnetic fields with arbitrary time dependence in a three-dimensional 
domain Q C M 3 . They read 

VxE(r,i) = -J^B(r,0, (la) 
VxH(r,t) = ^D(r,t) + J(r,t), (lb) 

with the spatial variable r e and the temporal variable t subject to bound- 
ary conditions specified at the domain boundary dfl and initial conditions 
specified at time to- The vectors of the electric field and flux density are 
denoted by E and D and the vectors of the magnetic field and flux density 
by H and B. The electric current density is denoted by J. However, we 
assume the domain to be source free and free of conductive currents (J = 0). 
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Furthermore, we assume heterogeneous, linear, isotropic, non-dispersive and 
time-independent materials in the constitutive relations 



The material parameters /i and e are the magnetic permeability and dielectric 
permittivity. At the domain boundary, we apply either electric (n x E = 0) 
or radiation boundary conditions (n x E = c/i(n x n x H)), where c denotes 
the local speed of light c = (e/i)" 1 / 2 . We also introduce the electromagnetic 
energy W contained in a volume V obtained by integrating the energy density 
w as 



This paper focuses on a general formulation of the DGM on non-regular 
hexahedral meshes as well as the projection of solutions during mesh adapta- 
tion. The issues of optimality of the projections and stability of the adaptive 
algorithm are addressed. Special emphasis is put on discussing the compu- 
tational efficiency To the best of our knowledge, this is the first publication 
dealing with dynamical hp-meshes for the Maxwell time-domain problem 
employing the DG method in three-dimensional space. 

As they are key aspects of adaptive and specifically /ip-adaptive methods, 
we will also address the issues of local error and smoothness estimation. 
This includes comments on the computational efficiency of the estimates. As 
estimators are not at the core of this article the discussion is, however, rather 
short. 

2. Discontinuous Galerkin discretization on non-regular hexahe- 
dral grids 

2. 1 . Discretization of space 

We perform a tesselation of the domain of interest Q into iV hexahedra 
Ti such that the tesselation T = Uili 71 is a polyhedral approximation of fi. 
The tesselation is not required to be regular, however, it is assumed to be 
derivable from a regular root tesselation T° by means of element bisections. 
The number of element bisections along each Cartesian coordinate, which is 
required to an obtain element i of T is referred to as the refinement levels 





(2a) 
(2b) 
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L x ,i, Lyj, L Zyi . As we allow for anisotropic bisecting the refinement levels of 
one element may differ. In case of isotropic refinement we simply use Lj. 
The intersection of two neighboring elements % D Tk is called their interface, 
which we denote as 2^. As we consider non- regular grids, every face Tj of a 
hexahedral element may be partitioned into several interfaces depending on 
the number of neighbors K such that Tj = Ufc=i-^ifc- The face orientation 
is described by the outward pointing unitary normal rij. The union of all 
faces is denoted as J 7 , and the internal faces J 7 \ dQ are denoted as J-" mt . 
Finally, the volume, area and length measures of elements, interfaces, faces 
and edges are referred to as \7l\, \Iik\, \J~j\ and \Td,i\, where d denotes any of 
the Cartesian coordinates. Every element of the tesselation T is related to a 
master element T = [— 1, l] 3 through the mapping Gi 

G, : f 7- : r H- r = (^kl + X . Q , + y . Q> + ^ , (4 ) 

where d i denotes the element center. 

2.2. General formulation 

Multiplying Maxwell's equations (CQ) by a test function ip(r) G if 1 (71), 
integrating over % and performing integration by parts yields 

j (V/i-^H- (W) x d 3 r + y ^(n x E) d 2 r = 0, (5a) 

j ^e^E+(V^)xHjd 3 r-| ^(nxH)d 2 r = 0, (5b) 
% d% 

where the explicit dependencies of r and t have been omitted. Equations (j5J) 
constitute the generic weak formulation of the time- dependent Maxwell's 
equations. In the following, we will replace the exact field solutions E and 
H by approximations using the discontinuous Galerkin framework. 

The space and time continuous electromagnetic field quantities are ap- 
proximated on T as 

N 

U(r,t)«U A (r,t)=0U i (r,t), (6) 
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where U G {E, H}. The element-local approximation Uj(r, t) reads 

U 4 (r,0=E u ?(*K(r) (7) 
v 

with the polynomial basis functions <p(r) and the time- dependent vector of 
coefficients 

nm = RM<M<M T i (8) 

representing the numerical degrees of freedom. The basis functions are de- 
fined with element-wise compact support, which is an essential property of 
DG methods 

(r) = P (r >' *l % (9) 
I 0, otherwise. 

We define the basis functions on the master element T and obtain the 
element-specific basis through the mapping Gi as 

<Pi = <P°Gr 1 (10) 

We employ Cartesian grids and tensor product basis functions of the form 



d S {a;,y,z} 

where p is a multi-index obtained from all pa = O..Pd- We denote by Pi = 
(P x ,i, Pyj, P Zj i) the local maximum approximation orders of element 71- The 
finite element space (FES) V p spanned by the basis functions is given by the 
tensor product of the respective one-dimensional spaces 

V p = ® g> Vf z with Vf d = span{^ d (r d ); < p d < P d } (12) 

The approximation may, thus, make use of different orders Pd in each of the 
coordinate directions, where we drop the subscript if they are equal. The 
basis functions are Legendre polynomials scaled such that Js|] 

Ti 10, otherwise. 

In the following the dependence of the spatial and temporal variable is not 
written down explicitly. 
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If now we were to substitute the exact electromagnetic field solution 
{E, H} for its approximation the surface integral term of (JS} cannot be evalu- 
ated straightforwardly at the internal faces J-" mt . This is due to the ambiguity 
of the DG approximation at any interface as a result of (j6]) and (Q. Weak 
continuity at internal faces is obtained locally by introducing numerical in- 
terface fluxes as 

J ip (n x IT) d 2 r, (14) 

BTi 

where U* is a unique interface value computed solely from Uj and Ufc , where 
7^ is a neighboring element. Common choices include centered and upwind 
fluxes. The centered interface value is given as 

U d,ik = 2 ( U d,k\Zik + U d,i\T ih ) ■ (15) 

Computing the upwind value is more involved. It is obtained as the exact 
solution of Maxwell's equations for piece-wise constant initial data after an 



infinitesimal time span, which is referred to as the Riemannian problem [17 
For the x-component of the electric and magnetic field at an interface with 
normal n 2 they read 

jri* (Yk\T ik Ex,k\T ik ~ Hy^\x tk ) + (Xi\X ik E X: i\x ik + H y ^\z ik ) fia\ 

h x,ik = v Tv ' ( lba v 

Ik\T lk + *i\T lk 

w * {Zk\x lk H Xtk \x lk + Ey tk \ Xlk ) + {Zi\i ik H Xii \ x . k - E y ^x ik ) . . 
H x,ik = y —7} • (lobj 

6k\T ik -T 6i\ Xik 

with the intrinsic impedance and admittance 

Other components are obtained by cycling the component indices and signs. 

Note that centered fluxes preserve the Hamiltonian structure of Maxwell's 
equations while this property does not carry over to the semi-discrete equa- 
tions when applying upwind fluxes due to the mixing of electric and magnetic 
quantities in f|T6|) . Consequently, an energy conservation property 0, Hj can 
be obtained with the centered flux formulation only, determining the kind 



of time integration schemes to be used as well [18[. Our implementation 
includes both flux types. 
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Having resolved the ambiguity at interfaces, we insert the approxima- 
tions into the weak formulation (jSJ) and follow the Galerkin procedure 
yielding the semi-discrete DG formulation 



x E*) d 2 r = 0, (18a) 



J Lf fi^H h - (V</f ) x E fc ) d 3 r + J <pf (n 
J (y* e^E, + (Vvf ) x H fc ) d 3 r - j <pf (n x H*,) d 2 r = 0, (18b) 



Ti dTi 



Vz = 1..N, Vgj = O..Pj. The volume integrals are referred to as the mass 
and stiffness terms, the surface integrals represent face fluxes. Note that no 
assumptions on the grid regularity have been made in the derivation. 

2.3. Employing non-regular grids containing high level hanging nodes 

Due to the strictly element-local support of the basis and test functions, 
the DGM is highly suited for the application on non-regular grids. The actual 
difference of the refinement levels Lj and of neighboring elements, i.e., the 
level of hanging nodes, plays a minor role as shown in the following. 

Inspecting equations (ITS]) it is seen that the mass and stiffness terms are 
not affected by the grid regularity as they are strictly local to the element %. 
The flux term, however, involves neighboring elements as well. Decomposing 
the surface integral into the six contributing face integrals 

J <pf (n x XJ* h ) d 2 r = itj <pf (n, x U*) d 2 r, (19) 
and considering centered fluxes for brevity each of these can be expressed as 



J cpf (n, x U,) d 2 r + J2 f ff K' >< U*) d 2 r 



X ik\j 



(20) 



Accounting for the kind of non-regular grids described above, i.e. grids ob- 
tained from a regular root tesselation, requires no more than summing up the 
contributions of all neighboring elements to the total flux. This is indepen- 
dent of the hanging node levels as well as the actual number of neighboring 
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elements. Inserting the approximation ([7]) into (120 j) yields 



1 
2 



x 



L Pi 



x u 



Pk 



(21) 
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Again, the first integral term does not depend on the grid regularity. As- 
suming rij to be aligned with the ^-coordinate and to point towards positive 
direction it amounts to 

j tfVfd 2 r = (22) 

due to the basis function scaling ([TBI . The second integral term can be 
expressed as 

J ipftft d 2 r = ^(l)^f(-l) | ^< fc dx | vZiVfikdy. (23) 

In this case the orthogonality property of the basis functions is lost due to 
non-identical supports of <fi and cpi-. We gather the terms (122]) and ( 123|) in 
the interior and exterior flux matrices F~ and F + . Following to the usual 
notation the sign indicates the evaluation from the interior and exterior side 
of the interface. Any non- regularity of the grid is now concealed within F + , 
which reduces to the standard form on regular grids. 

For high level hanging nodes the number of integrals to compute quickly 
becomes large, imposing a heavy computational burden if integration is per- 
formed at run time. However, as the integrals J d nd . ( p'di ( Pd d k^- r d m p3|) do 
not include the actual approximation but basis functions only, they can be 
precomputed analytically (making use of the master basis functions) and 
stored in tabulated form in the code. This has to be done for all combina- 
tions of pd and qd as well as for each possible edge overlap according to the 
respective difference in the refinement levels ALd (cf. Fig. [TJ). The number of 
possible overlaps grows as 2 ALd . We tabulated the integrals up to ALd — 6 
and for basis functions up to order six, yielding 247 matrices Ial of size 
7x7. In the isotropic refinement case ALd = 6 corresponds to one element 
interfacing with (2 6 ) 2 = 4096 neighbors. In the case of even larger differences 
in the refinement levels of neighboring elements, which are unlikely to occur 
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a numerical integration is invoked at run time. If the neighboring element 
has a smaller instead of higher refinement level the respective transposed 
matrix (Ial) t is applied. For upwind fluxes, the interior and exterior flux 
matrices do not change, however, they are applied to both, the electric and 
the magnetic field due to ffT6]) . 




Figure 1: Example of a non-matching interface. Black lines indicate edges of the root 
tesselation, gray lines indicate edges of refined elements. In the figure the interfaces are 
separated along the z-axis for a better visualization. The left hand root element has been 
refined several times, the right hand clement is at root level. The interface I connects an 
element of refinement levels (2, 3, L z ) with the root level element. The tick marks indicate 
possible locations for the imprint of elements of these refinement levels. The actual imprint 
on the root element face fills the first and sixth slab along the x- and y-axis, respectively. 
The interface II fills the respective second slab along the x-axis. 

In order to further enhance computational performance, all combinations 
of <p qd (±l)<p Pd (±l) and the integrals f^(-^ip gd )(p Pd df^ arising form the stiff- 
ness terms of (fT8|) are evaluated and tabulated as well. Precomputing the 
interface integrals maintains the high computational efficiency of the DG 
methods also for non-regular grids. Using matrix notation, the semi-discrete 
DG Maxwell equations f llSp read 



d_ (MM ( j(F~ + F+) Z -S + (F- + F+) \ fh 
dt\M € e) VS-(F- + F+) 7(F- + F+)/Z J \e 




where S and Z denotes the stiffness and impedance matrix. The matrix 
operator on the right hand side of (124)) represents a weak DG curl operator. 
Choosing 7 as either zero or one yields centered or upwind fluxes, respectively. 
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By applying centered fluxes the Hamiltonian structure of Maxwell's equations 
in continuum is preserved, whereas upwind fluxes lead to a mixed form. 
Symplectic explicit time integration can be applied in the former case but 



not in the latter one [18|. For examples of symplectic time integration for 
Maxwell's equations in the DG framework see, e.g., 0,0, In jaij upwind 



fluxes and Runge-Kutta schemes are applied for the time integration, where 
the latter one is concerned with Maxwell's equations. 

3. Local refinement techniques 

The adaptation techniques presented in the following are based on pro- 
jections between the finite element spaces introduced in (IT2"|) . The projection 



operators have been introduced in [16J, however, they are included for com 



pleteness. Also, we address the issue of stability in depth and amended this 
section with examples. 

The approximation to a given function / in the FES V p is obtained 
by performing an orthogonal projection. The projection is carried out in an 
element-wise manner, by means of the projection operator IP given by 

/, = E nP (/)^r = E#5r^ (25) 

p p WiiVitTi 

where (u, v)j: denotes the inner product J T uv dr on the element % with the 
associated 2-norm (u, u)j- = \\u\\ji. When applied successively to all elements 
and all components of given initial conditions of the electric field, E(t = to), 
and the magnetic field, H(t = to); the respective DG approximations and 
Hh are obtained. These approximations are optimal in the sense that the 
projection errors Sd = Ud — Udh are orthogonal to the space of basis functions 
V p 

(S d , tf) Ti = 0; Wp G [0, P], = <? o G~\ <<? G V p (26) 

3.1. h-Refinement 

As stated above /i-refinement is achieved by means of element bisections 
along the coordinate directions, where we allow for anisotropic refinements. 
The refined elements are referred to as the left and right hand side element 
7i and % with basis functions denoted as (p[ and <pl spanning the spaces V/" 
and V/ 2 in a full analogy to V p defined in f fT2|) . The approximation orders 
La and Rd in each child element do not have to be identical, neither are they 
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required to be equal to the respective order of the parent element. The 
direct sum of the spaces Vi and V r is denoted by V + 

V+ = Vi © V r . (27) 

In the following, the projection (1251) can be applied in order to project an 
approximation given in an element % to the FES associated with an /z-refmed 
or /i-reduced element. For h- refinement this yields 

(u,){ = IlftU^, (u^ = n;(U,) rr . (28) 

Due to the tensor product character of the basis, this can be expressed as 

= £X n i(tf)* = < * (vf)r. A tintfh, nj' (29) 

p p 

for the left and right child, respectively. If refinement is carried out along 
one coordinate only, e.g. x, this further simplifies to 

/x,l 

P:r Pa ' 

where 5 denotes the Kronecker delta. Note that above we loop over p x , 
whereas in ( 1291) the loop parameter is p. As, moreover, (tp*, $ x )% i vanishes 
for any p x < l x , we can limit the above loop to the range [l x ,Px,i], which 
reduces the number of addends to the minimum possible. 

For the merging of elements, the approximation within the parent ele- 
ment, Ti, is considered to be given piece- wise within its child elements. The 
projection reads 

uf = IPftU,), + (Ui) T ) Ti = TP{(Ui)i) Ti + ^((U,)^, (31) 
where the simplifications ( 129|) and ( 130]) apply. 
p-Refinement 

For the case of p-enrichments, the local FES are amended with the (P^+l) 
order basis functions 

Vf+^Vfu^f 1 }, (32) 

where any (non-zero) number of the local maximum approximation orders Pj 
may be increased. Also, an enrichment by more than one higher order basis 
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function is possible. Formally, we perform the orthogonal projection (I25p . 
however, due to the orthogonality property of the basis functions the coef- 
ficients u° " p remain unaltered under a projection from Vf to Vf +1 . Prac- 
tically, we simply extend the local vectors of coefficients Uj with the new 
coefficients uf +1 , which are initialized to zero. 

Conversely, for the case of a p-reduction, the local FES is reduced by 
discarding the P^-order basis functions 

Vf" 1 = Vf \ (33) 

Again, by virtue of the orthogonality, we find that the coefficients uf are 
deleted from the local vectors of coefficients while the coefficients u° " p ~ 1 
remain unaltered. 

We denote by IT7- the projection of the global approximation (E^, H/J 
from the current discretization to another one obtained by local h- and p- 
adaptations. 

3.3. Optimality, Efficiency and Stability 

3.3.1. Optimality 

An approximation with coefficients according to (125]) is optimal in the 
sense of (|26|) . The approximations within refined and merged elements with 
coefficients obtained through the orthogonal projections (128]) and (I3T]) are, 
hence, optimal in the same sense. 

If L d > P d and > Pd holds true for all d, the FES V is a subspace 
of V + (cf. (T27]) ) and every function of V is representable in V + but not vice 
versa. In this case, a given approximation is exactly represented within an 
element under /i-refinement but not under /i-reduction. See Fig. [2] for an 
example. 

3.3.2. Efficiency 

Since the projections (128]) for performing /i-refinement are independent 
of the actual approximation, we also tabulated the projection operators n| d 
and II£ d (expressed in master basis functions) yielding the matrix operators 
(IIi) + and (II r ) + , where the superscript denotes that the refinement level L 
is increased. Accordingly, we make use of the matrix operators (III) - and 
(n r ) _ for evaluating the projections of (|3T]) in the case of element merging. 
The matrix operators are related as 

(no- = 2(n+) T , (n r )- = 2(n r ) T . (34) 
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Figure 2: Projection based grid refinement and coarsening in one dimension. In (a) the 
projection of a given approximation (gray) to the left and right hand side child elements 
and the respective obtained approximations are shown (dashed, red/blue). If L x > P x 
and R x > P x , the approximations of the parent and child elements agree point- wise. The 
projection to a merged element shown in (b) can, in general, not be exact due to the 
discontinuity. 



This allows for the computation of the approximations within adapted ele- 
ments by means of efficient matrix-vector multiplications. As all projection 
matrices are triangular the evaluation can be carried out as an in-place op- 
eration requiring no allocation of temporary memory. 

3.3.3. Stability 

The global approximation associated with an adapted grid is computed as 
(Uj-Eh, Hj-Hh)- It can be considered as initial conditions applied on the new 
discretization obtained by performing the refinement operations. Assuming 
stability of the time stepping scheme (cf. (1, @, 0|), it is sufficient to show that 
the application of the projection operators at some time t* does not increase 
the electromagnetic energy associated with the approximate DG solution, 
i.e., 

W h (E h (f), H h (0) > W h (U r E h (t*), n r H h (t*)). (35) 

In this case it follows W h (t )) > W h {t*) > W h {T) and, thus, stability of the 
adaptive scheme. 

Following (J3J) the energy associated with element T% is given as 

Wi = ! \ (eE? + ^) d 3 r = U% (e,||e,||l + /^h,^) . (36) 
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As a consequence of (f2"9"|) . it is sufficient to show that the energy (135)1 is 
non-increasing during any adaptation involving one coordinate only. 

h- Refinement. For the following discussion of stability it is assumed that re- 
finement is carried out along the x- coordinate. Also we assume the maximum 
approximation orders L, R and P to be identical. It is clarified later, that 
this does not pose a restriction to the general validity of the results. 

In the case of /i-refinement, the operators 11+ and 11+ project from the 
space V to the larger space V + . Following the argument of paragraph 13.3.11 
on optimality, any function defined in the space V is exactly represented in 
V + . The conservation of the discrete energy is a direct consequence as the 
approximation in the parent and child elements are point-wise identicall. 

We find the following relation for the 2-norms of the respective local 
vectors of coefficients 

Wi = (Wi)i + (Wi) r 

^feNii + AiiMj!) = TfelK^illi + ^ll^illi + ^IK^rllS + AiilK^rlll) 

2(^11^111 + ^11^111) = e.dKeOilll + ||(ei) r ||i)+ ^(H (h^illl + HO^Hl). 

(37) 

The exemplary parent element approximation plotted in Fig. [2^ has a max- 
imum order of P = 6 with all coefficients equal to one. The coefficients 
of the child element approximations and the square values of their 2-norms 
are given in Tab. [TJ If the vector u is considered to be either the vector of 
coefficients of the electric field e or the magnetic field h the result agrees 
with (J37D. 





M 


Ml 


u 2 


u 3 


M4 


M 5 


M 6 


HIS 


U(x) 
U{x) x 
U{x) t 


1.0000 
0.2574 
1.7426 


1.0000 
-0.1355 
1.5631 


1.0000 
-0.2606 
1.6819 


1.0000 
-0.1446 
2.0399 


1.0000 
0.4276 
1.0494 


1.0000 
-0.1563 
0.2181 


1.0000 
0.0156 
0.0156 


7.0000 
0.3808 
13.6192 



Tabic 1: Parent and child element coefficients of the function plotted in Fig. [2^ 



The /i-coarsening operators Ilj - and 11" project a function from the space 
V + to the smaller space V. Since V is a subspace of V + , it is immediately 
concluded that, in general, energy is lost during the coarsening process. The 
discrete energy can only be preserved if the union of the left and right hand 



2 Identical material properties are assumed for the parent and child elements. 
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functions is an element of the smaller space V. Starting with the coefficients 
of the child elements, given in Tab. (TJ the parent element coefficients are 
exactly recovered and the discrete energy is preserved. 

In particular, it can be shown from algebraic properties of the projection 
matrices, that the discrete energy for arbitrary fine grid coefficients is always 
non-increasing during /i-coarsening. First, the (Pj x 2P^) projection matrix 
II - is defined as 

ir = (V nA (38) 



The coefficients of the child elements are gathered in one vector (uj) + 

K) + = ■ (39) 



u 



i It 



(41) 



Then, the coefficients of the parent element are given as 

u, = rr( Ui ) + , (40) 

which is equivalent to Eqn. ( l3~Tj) . Using (136]) and ( I57)) . the following must 
hold true in order to guarantee a non-increasing discrete energy 

2Klli k IIK)+||| 

2uJ U! k (K)+) T (u,)+ 

2 (n"K)+) T n-(u,)+ < (K)+) T (u,)+ 

(( Ul )+) T (n-) T n-(» i )+ ^ i 

(( Ul )+) T ( Ul )+ - 2- 

In order to fulfill this it is sufficient to demand 

max{ei g ((n-) T rr)} < 1 (42) 

Since the matrix ( (n _ ) T n~) has the P^-times degenerated eigenvalues 1/2 
and this is always fulfilled. 

The coefficients for the example shown in Fig. [2b are listed in Table EJ 
The sum of the 2-norms for the left and right hand child vectors of coefficients 
yields 1.9003 while twice the value obtained for the parent element evaluates 
to 1.8996. Thus, energy was lost during /i-coarsening. 
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u 


Ml 


U2 


^3 


Nil 


U(x) 
U(x)i 
U{x\ 


0.9500 
1.0000 
0.9000 


-0.0433 
0.2000 
-0.2000 


-0.2073 
-0.1000 
-0.0100 


0.0498 
-0.0100 
0.0100 


0.9498 
1.0501 
0.8502 



Table 2: Parent and child element coefficients of the function plotted in Fig. [2Jd 

p-Refinement. In order to show stability of the p-adaptation we again con- 
sider the energy stored in an element given by (1361) . In the case of p- 
enrichment, the local vectors of DoF are extended by the coefficients cor- 
responding to the (P + l)-order basis functions. Since these coefficients are 
initialized to zero, it holds true 

\\(n) P \\l = \\{uf\\l + = || (uo..p) p+1 Hi + ||(up +1 ) P+1 ||1 = \\(u) P+1 \\l (43) 

The discrete energy is exactly conserved. 

In the case of a p-reduction, the coefficients assigned to the highest order 
basis functions are removed from the vectors of DoF. Consequently, it holds 
true 

\\(u) p \\l = ||(uo..p-i) P ||^ + \\(u P ) p \\l < HK.P-ini! = IKu^-H, (44) 

and the discrete energy is either preserved or otherwise reduced. If the deci- 
sion for reducing the order is a correct one, the highest order coefficients are 
small and the induced energy loss is small. 

In the discussion of stability for /i-adaptations it was assumed that the 
maximum approximation orders L, R and P are identical. After showing that 
p-adaptation does not increase the electromagnetic energy either, it can be 
concluded that this assumption does not restrict the validity of the results 
obtained as the problem can be reduced to performing h- and p-adaptations 
sequentially. 

Finally, we will make some remarks on mixed h- and p-adaptations. If an 
/i-coarsening goes along with a p-enrichment, the latter should be performed 
first. It amounts to a projection to the superspace V p+1 . Projection errors are 
limited to the projection onto the space V^ 1 . If the projections are carried 
out in reversed order the final space is V^ +1 as well, however, the second 
projection from adds zeros to the local vector of DoF only, resulting in 
an increased overall projection error. On the contrary, if an /i-refinement 
goes along with a p-reduction, the former should be carried out first for the 
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same reason. However, the two adaptations can be performed simultaneously 
by employing projection matrices of the size {P^ cw x -P^ ld ). 

4. Automatic hp- adapt ivity: Error and smoothness estimation 

In the preceding section we presented efficient techniques for performing 
local h- and p- adapt at ions, which can be considered as the basic toolbox 
within the larger frame of an Zip-adaptive DG method. In this section we 
address the critical issues of locally estimating the approximation error and 
solution smoothness, which is required for driving the adaptation process. 
Note that as we are interested in estimating the smoothness of the actual 
solution, the discontinuous nature of the DG approximation is not a concern. 
In fact, we will show how the discontinuities at element boundaries can be 
exploited for constructing a smoothness indicator. 

In order to perform an adaptation of an hp-mesh two steps have to be 
carried out. First, the elements requiring adaptation have to be identified. 
This is achieved by an element-wise error estimation. If the estimate exceeds 
some tolerance the element requires refinement, elements having a small error 
are eligible for coarsening. In a second step, the kind of adaptation (h or 
p) suitable for the respective element has to be decided upon. Here, we 
distinguish between the refinement and the coarsening case. 

For the coarsening case the best option is obtained by consecutively test- 



ing possible de- refinements from a set of candidates |20|, |21[ . The set contains 
at least the candidates obtained by /i-reducing the refinement level to Lj — 1 
and by reducing the polynomial order to Pi — 1 but larger sets of candidates 
are possible as well. The list of candidates can be extended hierarchically. 
If, e.g., reducing the polynomial order yields an approximation still fulfilling 
the accuracy requirements, subsequent candidates can be generated as long 
as the accuracy demands are met. 



An extension to the refinement case is possible as shown in |2CH22| , how- 
ever, at the price of computing a globally h- and p-refined solution. We pur- 
sued a different approach based on a local smoothness indicator (cf. Sec. 14. 2p . 

4-1. Error Estimation 

Error estimation is addressed in a large number of publications out of 



which we refer to the introductions |23|. 24 and references therein. In the 



context of this paper, we focus on |25| , where a relation between the size of 
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the jumps of some quantity at element boundaries 



Phk = U AXik - U k]Iik (45) 

and the residual within in the respective element is derived for a model 
problem. Following this idea, we obtained a similar result for Maxwell's 
equations, which we apply as an error estimate. This is a valuable tool 
for the development of an adaptive DG method because the evaluation of 
the jumps is a computationally inexpensive operation. Moreover, it can be 
incorporated with the computation of the fluxes, which renders the extra 
costs negligible. The derivation is outlined in the following. 

The starting point is equation (I18ap . where we integrate the second term 
of the volume integral by parts to obtain 

J if> (JL^H h d 3 r + J i> (n x (E* - E h )) d 2 r + J ^(V x E,,)d 3 r = 0, (46) 

where we returned to denoting the test function by ip (cf. (jSJ)) for emphasizing 
the freedom of choice. Taking ip = 1 we obtain 



R, = J R h d r = J n x (E h - E* ) d 2 r, (47) 

with the residual 

R h = ^H h + (VxEfc), (48) 
and the element error estimate 

% = W\Tv (49) 

The global error estimate is obtained as £ — ^2i£i- By inserting either of 
the fluxes (TT5]1 or ( fT6l) . i.e. centered or upwind, in ( ]4Tjl we obtain 

Rf = ^Z)/ nx I E k d2r ' ( 50a ) 

R " P = \ E / n X d E l^ + Z Mnd d 2 r, (50b) 
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where k numbers all interfaces of element i. This establishes the link between 
the jump sizes and the residual. Note, that there is a consistency with the 
observation that in homogeneous regions the jumps disappear if the numerical 
solution is exact. 

Figure [3] shows results of tests we conducted for investigating the esti- 
mate performance. The tests were carried out in a one-dimensional domain 
using a Gaussian and a trapezoidal wave form as examples for a smooth so- 
lution and non-smooth solution. Bold lines in Fig. [3] correspond to the global 
approximation error in the L 2 -norm, dashed lines represent the global error 
estimate £. 

The estimate performs well for the Gaussian wave form. It tends to over- 
estimate the error, however, we consider the discrepancies to be acceptable, 
especially given that it reproduces the correct trend under grid refinement. 
The situation is not as good for the trapezoidal wave form, where the overes- 
timation of the error is more significant. Nevertheless, also for this example 
the trend under grid refinement is correct. 




Figure 3: Global approximation error for a Gaussian and a trapezoidal wave form under 
uniform grid refinement. Bold lines correspond to the L 2 -norm, dashed lines indicate the 
estimate. 



4-2. Smoothness Estimation 

Once an element is marked for refinement, it has to be decided upon the 
kind of adaptation to perform. This /ip-decision is based on the solution 
smoothness. If the approximation within the element under consideration 
is sufficiently smooth we perform p-refinement in order to obtain spectral 
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convergence [26|, |27j|. Otherwise, we choose h- refinement. Hence, the dis- 
cretization error and the achievable order of convergence critically depend 
on correct Zip-decisions. This in turn makes reliable smoothness estimators 
a necessity for adaptive codes. 

The subject of smoothness estimation has a good coverage in the litera- 



ture as well. In one of the first publications on adaptive DG methods (28 
the authors propose to estimate the local solution regularity based on the 
decay rate of the local error. As the local error is unknown and subject to 
estimation itself the adaptation procedure critically hinges on the error es- 



timate. More recent works such as (2j|-[34( attempt to estimate the solution 



smoothness based on a variety of properties of the local solution, but they 



do not include the error estimate. Sec 29] also for a more complete overview 



of different smoothness estimation concepts. 

A first family of popular methods estimates the local analyticity by pro- 
jecting the solution to a set of orthogonal functions (e.g. Legendre polyno- 
mials) and investigates the decay of the coefficients in on way or another. 



This approach is being followed by [29|, l3lH33l| , where the latter one builds 



on top of |3l| and obtains improved results especially when the sequence 
of coefficients exhibits a pronounced odd-even behavior. As we express the 
approximation as a Legendre series within each element the projection step 



would not be required. In [29j, |34j the authors attempt to estimate the local 



Sobolev regularity index directly, where the former one requires the repre- 
sentation of the solution as a Legendre series as well whereas the latter one is 
a novel idea based on continuous Sobolev embeddings and does not require 
a series expansion. Finally, smoothness indicators can also be built from su- 
per convergence properties of the DG method 30|,|35|, which will be described 
in more details below. 

The first family of smoothness estimates are mainly applied in the context 
of shock capturing, specifically for controlling artificial viscosity within high 
order DG simulations. In this context they proved to reliably achieve a good 
performance. In our experience, however, they are rather not suited for 
controlling an Zip-decision as they require a minimum number of coefficients 
for estimating the decay rate, which might not be available for low order 
elements. 

We adapted the smoothness indicator for hyperbolic conservation laws 



30 



which exploits the difference in the convergence rate of interface jumps for 
smooth and non-smooth solutions. From our experience this indicator works 
very robust and, most importantly, this remains valid down to very low 
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polynomial orders of one and even zero. This coincides with the experiences 



reported in |36|, |37(]. It should be noted, however, that no estimated value for 



the local analyticity or regularity is obtained. We plan to test the estimate 



341 ] since it appears to be working down to low orders while obtaining an 
estimated regularity index at the same time. For the time being, however, 
we consider an adapted version of the smoothness indicator 3^ ■ 



For a scalar quantity Q that can be a solution component or a derived 



quantity as well, it holds true [35 

/ (Q l -Q)d 2 r = 0(h 2p+1 ), (51) 



IX 



where X° fc ut denotes an outflow boundary regarding the quantity Q. Follow- 
ing 3(| the super convergence property (l5"Tj) can be exploited for constructing 



a smoothness indicator. To this end we consider 

/ [Q], fc d 2 r= f (Q l -Q k )d 2 r= f (Q l -Q)d 2 r + [ (Q-Q k )d 2 r. (52) 

As a result of (I5ip the last integral converges as 0(h 2(jp+1 ^) and the left hand 
side expression is 0(h p+2 ) in regions of smooth solution. If, however, Q is 
non-smooth in the vicinity of Ii k then both right hand side integrals are only 
O(h). Normalizing ( 152"j) by an average convergence rate and the L 2 -norm of 
the considered quantity on % yields a smoothness indicator 

|J^[Q,,Jd 2 r| 

hk ~ h(p+m\z lk \\\Q t \W (53) 

where h is a characteristic measure of the size of element i. We use the length 
of the edge normal to the interface X ik , which for the one- dimensional case 
is simply the element length. For Xj > 1 the solution on Xj is considered to 
be non-smooth and smooth otherwise with /i-refinement being performed in 
the former case and p-refinement in the latter one. As the indicator is ap- 
plied at element faces, it allows for indicating the solution smoothness along 
each coordinate separately which can be expoited for driving anisotropic hp- 
refinement. 

We use components of the Poynting vector S = E x H as the considered 
quantity Q. Then, inflow boundaries are recognized by S • < 0. The 
Poynting vector represents the energy flux density. As a consequence of this 



21 



particular choice, the nominator of ( 15 3 p is large for big jumps of the energy 
flux density across the element boundary. This is likely to indicate a low 
regularity of the local field solution and allows for interpreting the smoothness 
indicator based on the underlying physics. As this interpretation does not 
depend on the element order, we assume that this explains our observation 
of a good performance also for low orders (see Example 15.11) . 

5. Application examples 

In this section we will present two examples of the fully automatic hp- 
adaptive solution of wave propagation problems in one-dimensional space 
and one example which proves the capability of our implementation to han- 
dle large problems in three-dimensional space including thousands of mesh 
adaptations. However, in the latter example we drive the adaptation using 
a simple energy density criterion for the practical reason that the implemen- 
tation of the estimates for the three-dimensional case is not yet completed. 

5.1. Automatic adaptation in one- dimensional space 

The adaptation strategy described above is being tested with a Gaussian 
and a trapezoidal wave packet in a one- dimensional setting as shown in Fig. [4j 
The error tolerances are 10 -3 for the former and 10 _1 for the latter case. 
The gray dashed lines depict the position of the grid points and the red 
circles indicate the approximation order employed for the respective element 
divided by ten. For the Gaussian packet the adaptation algorithm chooses 
medium sized to big elements and the preset maximum approximation order 
of five in the vicinity of the packet. For the non-smooth trapezoidal packet, 
the maximum approximation order chosen by the algorithm throughout the 
simulation is two. In the vicinity of the pulse edges a high degree of h- 
refinement is applied, thus showing the desired behavior for the second packet 
as well. 

Figure |5] shows a comparision of the L 2 -error achieved for varying numbers 
of DoF using hp- adapt ivity (corresponding to different error tolerances) and 
fixed meshes of uniform polynomial order. Using hp-meshes the number of 
DoF required for obtaining a certain accuracy is clearly reduced. The achiev- 
able gain by using /ip-adaptivity, however, strongly depends on the problem 
under consideration. For the examples considered here, e.g., enlarging the 
domain while preserving the size of the wave packet will increase the sepa- 
ration between the curves corresponding to the adaptive and non-adaptive 
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solutions in favor of the Zip-solution and vice versa. In other words, the more 
pronounced the multiscale character of the problem is the higher is the gain 
by applying /ip-adaptivity. 

5.2. Proof of feasibility example in three-dimensional space 

In this section we consider the backscattering of a wide-band electro- 
magnetic pulse from a passive radar reflector. As stated above, the imple- 
mentation of the error and smoothness estimates for applications in three- 
dimensional space is underway, and we resorted to applying a physics based 
criterion for driving the mesh adaptation process. The purpose of this exam- 



ple is to demonstrate the ability of our implementation [38| to handle large 
scale problems using hp-meshes. This involves the handling of meshes with 
hanging nodes and thousands of adaptations on the element level, which are 
carried out by means of the efficient techniques presented above. 

In this example a radar reflector is illuminated off-center by a horn an- 
tenna, which emits a Gaussian-modulated sinusoidal pulse covering the fre- 
quency range from 20 to 30 GHz. The initial waveform is a TE 10 (transverse 
electric) wave. The setup is depicted in Fig. [6j where the antenna is shown 
in a cut view along with contour plots of the pulse at two instances in time 
(t = 1 ns and t = 10 ns). Table [3] lists the parameters and dimensions. We 
simulated the full scattering process starting from the excitation inside the 
waveguide to the recording of the reflected fields at the same position. The 
total propagation distance is about sixty wavelengths. 



Parameter Value 

Waveform mode TEio 
Waveform frequency range 20 - 30 GHz 

Waveguide type WR-62 

Waveguide width 15.8 mm 

Waveguide height 7.9 mm 

Horn width 39.7 mm 

Horn height 29.0 mm 

Horn depth 67.5 mm 

Table 3: Setup parameters and dimensions 

We chose a maximum /i-refinement level of two and the local element order 
to vary in between zero and four. The local energy density, introduced in ([3]), 
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serves as the criterion for controlling the adaptation procedure. Denoting by 
Wi the average energy density of element i and by Wi = Wi/w the normalized 
energy density with w = max {vJi}, we assigned the local refinement levels 

i 

according to u>j < 0.55 : L, = 0; G [0.55, 5) : Lj = 1; 5 < : Li — 2 and 
polynomial orders as < 0.55 : Pj = 0; Wj G [0.55,5) : p = G [5,25) : 
P = 2; ^ G [25, 35) : P = 3; w t G [35, 45) : P = 4 with 5 = 0.01. 

The initial discretization consisted of 45 x 35 x 100 = 157, 500 elements. 
During the simulation the number of elements varies and grew strongly after 
scattering from the reflector took place, when it reached close to 800,000 
elements corresponding to slightly more than 55 million DoF. For compar- 
ison, we note that employing the finest mesh resolution globally as well as 
fourth order approximations uniformly would lead to approximately 7.5 bil- 
lion (10 9 ) DoF. This corresponds to a factor of approximately 130 in terms 
of memory savings. We emphasize that the simulations were carried out 
on a single machine. The implementation takes full advantage of multi-core 
capabilities through OpenMP parallelization. The numerous run-time mem- 
ory allocations and deallocations are handled through a specialized memory 
management library based on memory blocking, which we implemented for 
supporting the main code 39 . 

Figure [7] depicts cut-views of the y-component of the electric field and 
the respective hp-mesh at three instances in time. Note that the scaling of 
the electric field differs for every time instance, which is necessary to allow 
for a visual inspection. The enlargement shows details of the computational 
grid. All elements are of hexahedral kind, however, we make use of the 
common tensor product visualization technique (cf. [40|,|4l|) using embedded 
tetrahedra for displaying the three tensor product orders (out of which only 
P x and P z are visible in the depicted x — z-plane). As we employed isotropic 
h- as well as p-refinement all tetrahedra associated with one element share 
the same color. Figure M shows plots of the outgoing and reflected waveform 
recorded along the waveguide center. The blue dashed line was obtained with 
the commercial CST Microwave Studio software 42j on a very fine mesh and 
serves as a cross comparison result. 



6. Conclusions 

We presented a discontinuous Galerkin formulation for non-regular hex- 
ahedral meshes and showed that hanging nodes of high level can easily be 
included into the framework. In fact, any non- regularity of the grid can 
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be included in a single term reflecting the contribution of neighboring ele- 
ments to the local interface flux. We demonstrated that the method can 
be implemented such that it maintains its computational efficiency also on 
non-regular and locally refined meshes as long as the mesh is derived from 
a regular root tesselation by means of element bisections. This is achieved 
by extensive tabulations of flux and projection matrices, which are obtained 
through (analytical) precomputations of integral terms. 

We also presented local refinement techniques for h- and p-refinements, 
which are based on projections between finite element spaces. These projec- 
tions were shown to guarantee minimal projection errors in the L 2 -sense and 
to lead to an overall stable time-domain scheme. 

Local error and smoothness estimates have been addressed, both of them 
relate to the size of the interface jumps of the DG solution. We considered 
the simulation of a smooth and a non-smooth waveform in a one-dimensional 
domain for validating the error and smoothness estimates. 

As an application example in three-dimensional space the backscattering 
of a broadband waveform from a radar reflector was considered. In this 
example the total wave propagation distance corresponds to approximately 
sixty wavelengths involving thousand of local mesh adaptations. As the 
implementation of the derived error and smoothness estimates for three- 
dimensional problems is subject of ongoing work, we chose to drive the grid 
adaptation using the energy density as refinement indicator. Crosschecking 
with a result obtained using a commercial software package showed good 
agreement. 
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Figure 4: Simulation of a Gaussian and trapezoidal wave packet using automated hp- 
adaptation. The electric and magnetic field arc plotted in blue and green, respectively. 
Gray dashed lines indicate grid node positions, red circles indicate the polynomial order 
of the respective element (divided by ten). The polynomial order is bound in between 
one and five. The initial, an intermediate and the final solutions are shown from top to 
bottom. For the Gaussian packet, the /i-refinement level does not exceed two, while taking 
full advantage of the highest order in the packet region. For the trapezoidal packet, an 
order of two is not exceeded. However, the algorithm makes use of /i-refinement levels up 
to six. 3Q 
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Figure 5: L 2 -crror of the solution vs. number of DoF for Gaussian (left) and trapezoidal 
wave packet (right) as depicted in Fig. [4] For the adaptive simulations the number of DoF 
corresponds to the mean value of all time steps. 




Figure 6: Scattering setup consisting of a horn antenna (in cut view) illuminating a radar 
reflector. The antenna emits a broadband electromagnetic waveform, whose electric field 
contours are depicted at two instances in time. 
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Figure 7: Visualizations of the y-component of the electric field (top panel) and the 
computational grid (bottom panel) at three instances in time. The enlargement shows 
details of the grid. Wc employ hcxahedral elements for the computation but make use 
of embedded tetrahedra for displaying the tensor product orders in the grid view. As 
isotropic p-refinement was employed in this examples all tetrahedra associated with one 
element share a common color. Note that different scalings are used for the time instances 
in the top panel. 
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